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Abstract 

We compute the spectral densities of T**" and J'^ in high temperature QCD plasmas at small frequency and 
momentum. uj,k ^ g'^T. The leading log Boltzmann equation is reformulated as a Fokker Planck equation 
with non-trivial boundary conditions, and the resulting partial differential equation is solved numerically 
in momentum space. The spectral densities of the current, shear, sound, and bulk channels exhibit a 
smooth transition from free streaming quasi-particles to ideal hydrodynamics. This transition is analyzed 
with conformal and non-conformal second order hydrodynamics, and a second order diffusion equation. 
We determine all of the second order transport coefheients which characterize the linear response in the 
hydrodynamic regime. 
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I. INTRODUCTION 



In relativistic heavy ion collisions a non-abelian plasma is formed which rapidly expands and 
evolves. There is a growing body of evidence from the Relativistic Heavy Ion Collider (RHIC) 
that this material equilibrates and can be characterized with viscous hydrodynamics [H, Indeed, 
viscous simulations of RHIC events indicate that the shear viscosity to entropy is remarkably small 



- < OAh , 
s 



(1.1) 



and the preferred value is close to the AdS/CFT prediction, rj/s = H/Att 10j,llll- From a theoretical 
perspective it is important to corroborate this phenomenological conclusion from first principle 
lattice QCD simulations. It is also important for lattice simulations to characterize what this 
nearly perfect fluid is like. For instance, at weak coupling there is a clear distinction between the 
inverse temperature ~ 1/T and the typical relaxation time ~ l/^^r 13, lid]) and consequently 
kinetic theory can be used to describe the real time dynamics of the plasma. In contrast, at least 
for the strongly coupled gauge theories which can be studied with AdS/CFT, there is no distinction 
between these time scales, and a quasi-particle description is hopeless Il5l | . This is reflected by 
the fact that there is no visible transport peak in strongly coupled spectral densities 0, [13] • 

Lattice QCD simulations at finite temperature can measure Euclidean correlators of conserved 
currents, ( J(t, £c) J(0, 0)) . These correlators are related to the spectral density of the corresponding 
operators by an integral transform 1S| 



d3cce^*'-^(J(r,cc)J(0,0)) 



^JJ(, , cosh(a;(r - 1/2T)) 
2^^ sinh(^/2r) 



(1.2) 



Generally, only the gross features of the spectral density can be determined from Euclidean mea- 
surements 



ig 



20|. Nevertheless, the urgent need for non-perturbative information about transport 



and other real time processes is evident from the strenuous analysis of available lattice data by 



several groups 2114261]. In an effort to characterize the response more completely, a number of 



lattice groups have begun to simulate the current-current correlators at finite spatial momentum 
271 . |28| . which reveals the real time information hidden in the Euclidean data as far as possible. 
The primary goal of this work is to calculate the spectral densities at finite uj and k at weak cou- 
pling for all possible combinations of T^*^ and J^. At zero k the weakly coupled shear, bulk, and 
current spectral weights have been determined previously in a full leading order calculation 29j, |30|. 
The spectral weights computed in this work exhibit in detail the transition from free-streaming 
quasi-particles to hydrodynamics. Ultimately, these perturbative spectral densities can be com- 
pared to the AdS/CFT and to the lattice, although a fair comparison to the lattice data requires 
a model for the high frequency continuum [311]. 



At weak coupling the kinetics of hot plasmas consists of several processes 13, ISJ]. First, 
the streaming of hard particles creates a random soft background field which causes momentum 
diffusion of the instigating hard particles. Second, there are collisions between the hard particles 
which cause an 0(1) change in a particle's direction. Finally, there is collinear bremsstrahlung 
which plays a particularly important role in the equilibration of the highest momentum modes 
331 ]. In the current work we will reanalyze the Boltzmann equation and make the diffusive process 
more explicit by reformulating the evolution as a Fokker-Planck equation. This analysis is limited 
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to a leading log approximation where the hard collisions and collinear bremsstrahlung are neglected 
for typical particles^. 

The QCD kinetic theory described above has been used to compute the transport coefficients in 



high temperature gauge theories 344361 1 , and to provide an intriguing, but incomplete, picture of 



equilibration in heavy ion collisions [3^, l37H39l| . Ultimately, one would like to simulate heavy ion 
event using the QCD Boltzmann equation, quantifying equilibration in the weakly coupled limit. 
While there has been considerable progress in understanding the dynamics of the soft background 



gauge fields out of equilibrium [4CH42I] , there have been only nascent steps in simulating the coupled 
particle field problem 43|, |4J]. As a by-product of the spectral weight calculation, we have made 
some (limited) progress in simulating the QCD kinetic theory. Traditionally, the shear viscosity 
is determined by minimizing a variational ansatz. However, when computing the spectral weights 
at finite fc, the variational ansatz would have to be two dimensional, and would have to capture 
a complicated structure reflecting the transition from Landau damping to hydrodynamics. In the 
two dimensional case it is easier to discretize momentum space and to solve the Fokker Planck 
equation directly. Finding the correct solution requires understanding the appropriate boundary 
conditions. In particular, there is an absorptive boundary condition at p = which accounts for 
the flux of particles from the temperature scale to the Debye scale and rapid number changing 



processes in the p — )• limit [361 ]. Section ITTl reviews this boundary condition and shows how gluon 
number equilibrates while conserving energy and momentum. Now that the problem of determining 
spectral densities has been reformulated as a definite initial value problem, the resulting numerical 
procedure can be used to simulate jet medium interactions in detail. In fact, it was this goal that 
led to the present work. 

Throughout, we will denote 4-vectors with capital letters P, Q and use p, q for their 3-vector 
components, Ep,Eq for their energy components, and p,q for |p|,|q|. Spatial 4-vectors and 3- 
vectors follow an analogous convention, e.g. X,Y and x,y, respectively. Our metric convention 
is [-,-|-,-|-,-|-], so that u^u^ = —1. We will notate the Bose and Fermi equilibrium distribution 
functions with = l/{e^^'^ — 1) and rip = YjiePl'^ A- 1), but will drop the BjF label when 
the appropriate statistics are clear from context. Momentum space integrals are abbreviated, 
4^/dV(2vr)3. 

II. LINEARIZED BOLTZMANN EQUATION 

Our starting point is the Boltzmann equation 

{dt-rv^-d^)f(t,x,p) = Q\Uv\. (2.1) 
which we will linearize around the equilibrium distribution with constant temperature To 

f{t,x,p) = np + 6f{t,x,p) , Up = ^^^^^ _ ^ . (2.2) 

The background is characterized by the energy density Cq, the pressure Vo, a squared sound speed 
c^, and the specific heat Tci,. Initially we will consider only pure glue theory and subsequently 



^ As discussed in the text, this is not exactly true. Bremmsstrhalung will determine the boundary conditions of the 
Fokker Plank operator as p — >■ 0, but otherwise can be neglected. 
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extend the analysis to include quarks in Section |Vl In a leading-log approximation, the only 
diagram is t— channel gluon exchange. Appendix |A] re- interprets the linearized Boltzmann equation 
(at leading log) as a Fokker-Planck equation, which is useful in subsequent analysis. Appendix [Xl 
extends the analysis of Refs. 



34, 



461 ] to general partial waves, and gives explicit expressions 



for the leading log gain terms. The gain terms do not affect the shear viscosity, since they only 
contribute to the i = and i = 1 partial waves. However, the gain terms do affect the spectral 
weights at finite k and oj which is the primary focus of the current work. 
Following Appendix [XJ the linearized Boltzmann equation can be written 



d 



{dt + Vp ■ dy,)6f = TfiA^ np{l + rip)—: 



d 



dp'- 



np{l + Up) 



+ gain terms , 



(2.3) 



where «4 is the drag coefficient of a high momentum gluon in the leading log approximation 
scheme [471, \^ 



dp 

d^ = -^^^' 



fJ-A = z log I 



Svr 



and the Debye mass for a pure glue theory is^ 

-2^- r np{l + np) _g^T^ 



2 9'Ca 
rujj = — : — u„ 



dA 



mo J 



N,. 



(2.4) 



(2.5) 



Without the gain terms, Eq. (|2.3p is a Fokker-Planck equation describing a random walk of the 
hard particles. 

In the diffusion process, the momentum-space current is given by 



jp = -T/iAnpil + Up) 



d_ 
dp 



np{l + Up) 



(2.6) 



and the work on the particles per time, per Degree Of Freedom (DOF), per volume can be found 
by multiplying both sides of Eq. (12. 3p by Ep and integrating over phase space 



dE 
d^ 



d'^p 
(2^ 



P-3p 



The momentum transfer (per time, per DOF, per volume) is similarly 



dP 

~dt 



d'^p . 



(2.7) 



(2.8) 



With these definitions and Appendix [XJ the gain terms are 



1 

^ain terms = — 
c,B 



1 (9 2 N 



dE 1 
'dt^^ 



d_ 
dp 



— np(l + Up) 



dP 

~dt 



(2. 



^ We follow an almost standard notation. The dimensions of the adjoint and fundamental representations are 
dA = — 1 and dp = A^c The Casimirs of the adjoint and fundamental are Ca ~ Nc and Cf = {N^ — l)/2Nc. 
The trace normalization of the adjoint and fundamental are Ta = Nc and Tf = 1/2. Ug = 2dA and i/q = 2c?f 
count the spin and color degrees of freedom for gluons and quarks respectively. 
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where for subsequent use we have defined 

" y (^"^^^ + = ^ , Cf^J -n,) = -. (2.10) 

In this form it is straightforward to verify that energy and momentum are conserved during the 
forward evolution of the hnearized Boltzmann equation described by Eq. (j2.3p and Eq. (j2.9p . What 
is less obvious is that gluon number is not conserved, as wih be discussed in the next section. 



A. Boundary conditions and particle flux to low momentum 



The resulting integral differential equation is ill-posed without boundary conditions. We will 
discuss the boundary conditions at low and high momentum respectively. 



1. Boundary conditions at zero momentum and number non-conservation 



To discuss the appropriate boundary conditions for the gluon number, consider the excess of 
soft gluons within a small ball of radius Ap ~ gT centered at p = 



/.p=Ap ^3 



fp=Ap rp2 
'p=0 

where here and below we define 



Sf{p)=np{l + np)x{p). (2.12) 

Since it is easy to emit a soft gluon it is easy to intuit that the appropriate boundary condition is 

X(p)lp^0 = 0. (2.13) 

As discussed in the introduction, Bremmstrahlung can be neglected for momenta of order ~ T in 
a leading log approximation. However, since the Breammstrahlung rate increases as the momentum 
is lowered, there is a momentum of order ~ gT where inelastic processes are important for any 
arbitrarily small coupling constant. Ref. [361 ] (see section E) estimates that the total rate for the 
hard particles to absorb (or emit) a gluon from the ball of radius Ap ~ gT is 

T\ti ~ g'T ^fip) ~ /T^ % . (2.14) 
Jo P Jo P 

The infrared divergence of the integral is cut off by the thermal mass of the radiated gluon m ~ gT, 
so that the total rate is of order g^T 



m 



rf^ai^y /T. (2.15) 



We are interested in the evolution of the system on a time scale, l/g'^Tlog(l/5f), which is large 
compared to the inverse radiation rate, 1/g^T. As we observe the evolution on this longer time scale, 
the soft Bremsstrahlung rate will rapidly maintain the equilbrium phase space distribution of soft 
~ gT gluons up to corrections of order of the ratio of these two time scales, ~ g. Thus, the excess 
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from equilibrium at small momentum, x{^p)^P-> is small at small momentum, xid^) = 0{g). 
The boundary condtion in Eq. (j2.13p is the leading result in the weak coupling limit g ^ 0. 
Even in a full leading order anaysis (where thermal masses are neglected on external lines), the 
boundary condition must still be adopted to avoid a divergent soft Bremmstrhalung rate [3Q]. A 
consequence of this boundary condition is that gluon number is not conserved during the Fokker 
Planck evolution as will be discussed in Section IIIBi 



2. Boundary conditions at high momentum 

At high momentum, the second order differential equation typically consists of an exponentially 
growing solution and an additional solution behaving at most as a polynomial. Clearly we should 
select the second solution as physically acceptable. 

To determine how to select the appropriate boundary conditions at high momentum, we reex- 
amine the Fokker Planck equation 

{dt + Vp-d^)6f{t,x,p) = -fiA{l + 2np)p-^ + TfiA^ldf. (2.16) 

The motion of the particle excess can be described alternatively by Langevin evolution with drag 
and random kicks ^(t) 

^ = -^^A{l + 27ip)f + e{t) , {em' it')) = 2TfiA6'^{t - t') . (2.17) 

At high momentum, we neglect the noise term Vp6f, set (1 + 2np) ~ 1, and keep only the drag 
term, yielding 

{dt + Vp • dx)6f{t,x,p) = -fiAP- ^ • (2.18) 

The resulting differential equation is now first order in derivatives and can be used to choose a 
boundary condition. Specifically we will discretize the first derivative 



P 



d6fit,x,p) _ 1 
dp p, 



,2 



) ) Pn ) ) </*p ) 

Ap 



(2.19) 



and use Eq. (j2.18p to solve for the first point off the momentum space grid, f{t,x,pn+i,9p,(pp). 
This selects the appropriate non-exponentially growing solution. 



B. Evolution of simple initial conditions and the flux of gluons 

If df{t,x,p) is given some initial condition, the resulting Fokker-Planck evolution will equi- 
librate, and 5f(t,x,p) will ultimately reach a form described by linearized hydrodynamics. 
Specifically, the system will be described by a temperature excess 6T(t, x) and flow velocity 
= {l,u'^{t,x)) that obey the equations of linearized hydrodynamics. The distribution func- 
tion will approach 

/eq(t, X,p) = ^^p.uit,^y^T^+sTit,a,)) _ l = + + i^ ^^xf^ + Yo^'^^^ '^^) ' ^^"^^^ 
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i.e. xiti^iP) will approach an equilibrium value 

Xec,uil{t,X,p) = p^^^^fi + ^u\t,x) . (2.21) 

It is instructive to examine how the total number of gluons, SNpp = J^5f{p), equilibrates 
during the Fokker-Planck evolution. Integrating both sides of Eq. (j2.3p yields two terms 



where ilp is an outward directed solid angle. The first term represents a diffusion flux of gluons 
with momenta of order p ~ T to momenta of order p ~ gT. The second term is the number of 
gluons disturbed from equilibrium per unit time by the random walk of the excess, 6f. In a given 
time step, this number disturbed is proportional to minus the work done on the excess by the 
bath, i.e. —dE/dt is the work done on the bath by the excess. In general, these two rates are 
different and number changes accordingly. In equilibrium, however, the two rates are equal and the 
excess number of gluons remains constant. Specifically, substituting the equilibrium distribution 
Eq. (|2.2ip . we find a net loss of excess gluons to the Debye scale, and a net gain of excess gluons 
due to the work on the bath 

Further information about the flux of particles from the temperature to the Debye scale is 
obtained by analyzing the small momentum limit of Eq. (j2.3p . (For simplicity, we will ignore any 
spatial dependence of 6f.) The expansion of x{p) near p ~ is characterized by its £ = and 
1=1 spherical harmonic components 



x{p) ^ 



dp 



=0 dp 



•p, (2.24) 

p=0 



where Xi is the Cartesian translation of Xim(p) in a spherical harmonic expansion. Substituting 
this form into Eq. (j2.3p , we notice a divergent rate as p — )■ 



^ , , 2Tp / dxi(t) 



p \ dp 



p=0 




(2.25) 



Because of this divergent rate, the slope at the origin will rapidly adjust itself to maintain the 
balance condition 



dXi{t) 



dp 



p=0 TflA^B dt 



1 dP , , 

(2.26) 



Thus the angular dependence of the flux (i.e. p^jp as p — t- 0) is determined by the momentum 
transfer to the hard particles by the bath. It is straightforward to verify that the equilibrium 
solution Eq. (j2.2ip satisfies this balance condition. 

Clearly the evolution of the excess is complicated, and the structure of Eq. (j2.3p is quite singular 
in the infrared. To gain some intuition about the solutions to this equation, we will show how an 
initial out of equilibrium distribution approaches equilibrium. For simplicity, we will consider two 
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out of equilibrium initial conditions. The first is spherically symmetric and independent of spatial 
coordinates, while the second initial condition is also spatially independent but is proportional to 
the £ = 1 spherical harmonic. Specifically, for the two cases we have 



x{t,p)=Xoo{p,t)Hoo{p) (easel), (2.27) 
x{t,p)=Xio{p,t)Hw{p) (case 2), (2.28) 

where Hqq{p) and Hiq(p) are the i = and £ = 1 spherical harmonics. For the initial condition at 
time to = 0, we take 

p^ripil + np)xoo{p, to) or p^Upil + np)xio(p, to)] oc ^ se-^P-^^'o)'/^.^ ^ ^2.29) 

with pq = 3To and o"^ = Tq. The numerical procedure to solve Eq. (j2.3p for the £ = and £ = 1 
partial waves is elaborated in detail in Appendix [Bj Briefly, Eq. (j2.3p (without the gain terms) is 
a parabolic differential equation. The momentum space is discretized, an implicit scheme is used 
to perform the update step, and the conjugate gradient algorithm is used to perform the matrix 
inversion. Fig. [U^a) and (b) show how the two initial conditions evolve as a function of time. At 
late times, the two initial conditions approach the equilibrium distribution Eq. (I2.2ip where the 
parameters 5T and n' are determined by the total energy and momentum in the initial state, e.g. 

" = ^;tv'o I (^ P + ■ (2.30) 

We have verified that Eq. (j2.22p and Eq. (j2.26p are satisfied during the evolution. Although the 
structure of the evolution equations is quite singular in the infrared, due to the boundary conditions 
(Eq. (j2.13p ) and the momentum balance condition (Eq. (j2.26p ). the final solutions are smooth and 
generally regular functions of momentum. 



III. SPECTRAL DENSITIES OF T^" 
A. Preliminaries 

Our goal is to compute all spectral functions of T^'^ as a function of u and k. To this end, we 
will compute the retarded Green functions 

/oo />oo . |- -, , / 

dt da;e+*'^*-*'^-^^(t)( r^'^(t,a;),r"'^(0,0) ) , (3.1) 
-oo J — oo 

and the associated spectral functions 

p^'•'"'^{uJ,k) = -2ImG^^°'^(w,fc) . (3.2) 



Here the ' indicates that the average is over the partition function in flat space [49[ . The distinction 
is necessary since the easiest procedure to calculate such correlators in kinetic theory is to disturb 
the theory with a weak external gravitational field. 

In the presence of a weak gravitational field, the stress tensor of the fluid in perfect equilibrium 

is 

TtqiX) = (e(r) + V{T))u^'{X)u''{X) + V{T)g^''{X) , (3.3) 



8 




FIG. 1; Evolution of an initial condition towards equilibrium in the linearized Boltzmann equation, (a) 
A spherically symmetric {(. = 0) initial condition approaching equilibrium at various time steps. Steps are 
in units of T/^a with = g^CA'rnjj\og{T/'mD)/STr. (b) An initial condition proportional to the first 
spherical harmonic Hiq{p) evolving towards equilibrium. In each plot the dotted lines show time steps of 
Q.2T/ jiA- The solid lines show times steps in units of 1.0r//i^. 



where = {1/ \J —gQQ{X),Q). Corrections to this equilibrium form can be found from the modifi- 
cations to the density matrix which is perturbed by the gravitational field. To linear order in the 
gravitational field, the action describing material in curved space time is 



.-int 



and the interaction Hamiltonian is 

Hint{t) = - J d'xC 
Then following the usual discussion of linear response, we have 

{T'"'{X))j,^^=T>^^'{X)-^ I d'Ye{X^-Y^)([T^'^{X),T"^{Y) 
In Fourier space we have for K = (oj, k) ^ 0, 



a0 



dh 



a/3 



(w, k) - ]-G^j^'^'^{u, k) Kp{uj, k) . 



(3.4) 



(3.5) 



(3.6) 



(3.7) 



h=0 



Thus to determine the spectral functions, we will turn on a weak gravitational field in the kinetic 
theory and determine the attending change in the distribution function and the stress tensor. 

To classify the relevant correlators, we will choose k along the z axis. Then the correlators can 
be classified according to their transformation properties under rotations around the z axis. In the 
next sections we will compute the following complete set of response functions: 



1. G^j^^^{u}, k) - Shear mode 
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2. G^j^^^{uj,k) - Sound mode 

3. G^^^{oj, k) - Tensor mode 

4. rifj^yTjapG^^'^^ {ijj , k) - Bulk mode 



B. Shear mode 



In this section we will compute the retarded correlator G^^^{uj, k) by turning on a gravitational 
field hzxii,z). Then we will analyze this correlator in the free-streaming and hydrodynamic ap- 
proximations to determine the high and low frequency limits respectively. The procedure (which is 
reasonably standard , 0, [s^]) will be described in detail in this section, and the same procedure 



will be used subsequently without explanation. 

In the presence of a gravitational field, the streaming term of the Boltzmann equation becomes 



5l|,l52l 



(^'sl^ - ^^^'^^aPx) = C[fM , (3.8) 

where 

= ^g^' {d,gp, + d,gf,p - dpgp,) . (3.9) 

In this equation df/dP^ should be understood as zero. Except in the bulk channel, all temporal 
components of the metric perturbation are zero 

9ij - Vij + hij . (3.10) 

Then we linearize the Boltzmann equation around the equilibrium distribution, which also depends 
on the background metric 

f{t,x,p) = n'l + 6f{t,x,p), nh = ^ = . (3.11) 

Substituting Eq. (j3.1ip into Eq. (j3.8p and linearizing with respect to 6f and hij leads to the 
following equation for 6f 

{dt + vp-d^)6f + np{l + np)^^dthij = C[5/,p] . (3.12) 

In Fourier space, this equation reads 

p^p^ d ( f^x(p)\ 

+ gain terms . (3.13) 

Without the gain terms, this is a linear elliptic partial differential equation which results from 
the steady state limit parabolic differential equation, as is usually the case. Traditionally, such 
differential equations are solved by matrix inversion, e.g. by the conjugate gradient method. We 
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will adopt this approach; in Appendix [Bl we introduce a spherical harmonic basis and discretize 
the radial momenta. Then the Fokker-Planck operator is realized with a straightforward second 
order difference scheme. This procedure, with only h^x 7^ for the shear channel, determines 
5f{uj,k)/hzx{'^,k) which can be used in subsequent analysis. 

After solving for 6f{uj, k), the energy- momentum tensor can be computed from kinetic theory 



rj-}Z 



(n^p + 5f{u,k)^ , (3.14) 



(27r)3 Ep 

^3p pZpX k) 



-Vohzxi'^^k) + 



(27r)3 Ep h;,x{(^,k) 



hzxiu},k). (3.15) 



Comparison with Eq. (j3.7|) shows that the term in square brackets is minus the retarded correla- 
tor, —G^^^^{uj,k). Fig. WlySi) shows the associated spectral density in the shear channel, and this 
correlator will be analyzed at low and high frequency in the remainder of this section. 

At low frequencies, hydrodynamics in an external gravitational field describes the behavior 
of this correlator. Hydrodynamics consists of the conservation laws together with a constituent 
relation 

V^T^^^ = , T^"' = T;^;;^! + TT^^ + n A^^ , (3.16) 

where T-f^^^ = {e{T)+V{T))u'^u'' +V{T)g'"' , A'^'" = gt^'+u'^u'' is the projector, and vr^'' is traceless. 
The strains n'^'^ and 11 are expanded in gradients, vr'^'" = ir^'^ + 112'^ + . . . and 11 = Hi + 112 + • • - j 
with 

vrf'^ = -rjaf"" , Ui = -QV^u^ . (3.17) 

Here denotes the covariant derivative, the brackets denote the symmetric, traceless, and spatial 
component of the bracketed tensor 

(V) = Ia^-A"^^ (^A^^ + A(s^ - V) ' (3-18) 

and a^"^ = 2(V^n'^). The second order theory which describes irt^'^ and H2 will be discussed 
in Section IIVI The weak gravitational field disturbs the energy momentum tensor away from 
equilibrium 

e{t,x) c^eo + eit,z), n'' ~ (1, n^(t, z)) , (3.19) 

where e and it* are first order in the perturbation. To first order in the field and disturbance the 
constituent relation becomes 

Tij = {6'^ - hij) + cl e5'^ - 2r? (d'u^) - CS'^dm' - r,dt {K,) - ^-^S'^dth , (3.20) 

where h = hn/S. The linearized equations of motion are 

dte + {eo + Vo)diu' = -^{eo + Vo)dth, (3.21) 
(eo + Vo)dtu' + djT^' = -Vo djhji . (3.22) 
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For the shear channel we specialize perturbation given above, with only hzx{t,z) non-zero. Then 
the equations of motion are easily solved in Fourier space 

e{uj,k) = 0, u'{uj,k) = 0, {co + Vo)u'' {oj , k) = ^'^'^ . (3.23) 

Substituting these results into stress tensor Eq. ()3.20p . we determine the retarded Green function 
in a first order hydrodynamic approximation 

2 

T'^oo, k) = -VoK^iu, k) - Cr^iuj, k)hU^, k) , Gl^'^iuj, k) = ""^"^ . (3.24) 

The imaginary part of this retarded Green function G^^^{ljJ, k) 



(cj and k small) , (3.25) 



describes the rapid behavior at small k and small uj seen in Fig. [2|^a). 
For high frequency, the free streaming Boltzmann equation, 

{dt + Vp ■ d^) 6f + np(l + np)^dthij = , (3.26) 

determines the spectral function. Again specializing the discussion to the shear mode where only 
hzxi^-, k) 7^ 0, we solve for 5f{oj, k). Then the stress tensor is 



-Vohzx{^^,k) + 



p Ep hzx{u},k) 



hzx{u},k), Sf{uj,k) 



pZpX -cjhzxnpjl + Up) 

Ep UJ — Vp ■ k + ie 



(3.27) 

Again the quantity in square brackets is the free streaming prediction for the response function 
—G^^^{uj,k), and we have introduced the +ie in order to specify the retarded response. Taking 
the imaginary part of the response function (using Im [1/(2; + ie)] = — 7r5(a;)), we determine the 
spectral density from the free theory 

(., k) .,vr3 -'f,_^) _ ^ k large) . (3.28) 



2uj 30 P V 

The free solution is shown as a dashed line in Fig. ^a) and agrees with the full result at large u 
and k except near the light cone. 

C. Sound mode 

In the sound channel, the only non-zero metric perturbation is /i^z(w,/c). The hydrodynamic 
prediction for 7^ is 

1 2 _ - n 3 

T'%u, k) = -Vohzziio, k) - -G'^'^iu, k)hzz{u;, k) , G^^^^(a;, k) = (e^ + Vo) 



2 5 ) ) ) ° ° _ ^2^2 _|_ iYgk'^oj ' 

(3.29) 

where = (|r/ -|- O/i^o + 'Po)- The free prediction is 
and is shown in Fig. [2|^b). 
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D. Tensor mode 



In the tensor channel, the only non-zero metric perturbation is hxy{io,k). The hydrodynamic 
prediction is then 

T^y{LO, k) = -Vo Kyiio, k) - GY^i^, k)Ky{Lo, k) , k) = -iurj , (3.31) 

while the free prediction is 



2oj ~ 120 k 



(3.32) 



E. Bulk mode 



To calculate the bulk spectral weight, considerably more care is required. In this case, the 
system is perturbed by a metric tensor of the following form 

g^,{X) = {l + H{X))T]^,. (3.33) 

If the gravitational perturbation is independent of time, hydrostatic equilibrium will be reached 
when T{x)^y —gQo{x) reaches a constant. Motivated by this observation, we show that for a 
conformally invariant theory, the distribution 

np{t,x,p) = j^^^j^^^^^^ ^ ^ (3.34) 

is a solution to the Boltzmann equation in the presence of the time dependent perturbation, where 

TH{X),/-goo{X) = Const , Tnix) = ^ (^1 - ^^^(^)) , (3-35) 

and Ufj{X) = {l/y^—gQQ{X), 0, 0, 0). This is true even if the temporal and spatial variations of 
H{X) are short compared to the mean free path, ~ 1/g'^T. Note that in the distribution function 
np{t, x,p), the combination P ■ U is 

- P ■ U{X) = 7^°^^!^ = Jp^{5,j + H{X)5,j)pi + m^{TH{X)) , (3.36) 

since P^'P''g^u = -m^{TH{X)). 

In the presence of a gravitational field and a nontrivial dispersion relation, the Boltzmann 
equation is 

- { - i^!^-^ - r^P^P'^^) fit, x,p) = C[f,p] . (3.37) 



The mass term might be more recognizable as a force term 

1 dm^{X) df _ dEpdf 
~ 2Ep dXt^ dPf, ~ ~~dxdp 



(3.38) 
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The need for this term when computing the bulk viscosity has been emphasized previously 531 ] . 
The mass depends on the distribution function, which in turn depends on space-time 32 1 



m\X)=g\(t>)CA4>{X). 



^(P°)27r<^(-p2)/p. 



(3.39) 



In equilibrium, (j) = . It is important to emphasize that in order for kinetic theory to be 

valid, the scale of variation of g^y{X) must be long compared to ~ ^/g^T which is the time it takes 
to establish the quasi-particle mass. Thus the gravitational field may be considered constant in 
Eq. 

We will substitute a trial solution 



f{t,x,p) 



npit,x,p)+5f{t,x,p 



(3.40) 

^ Therefore, 6f may be neglected 



and subsequently verify that 6f is of order (c^ — 1/3) ~ g 
when determining the gluon mass to leading order. The mass is then simply the time dependent 
equilibrium mass 



m\TH{X))^m\To) 



QT2 



H{X), 



(3.41) 



To 



where for pure glue, m'^{T) = g'^{T)CAT'^ /Q [32]. 

With this observation, we substitute Eq. (j3.40p into Eq. (|3.37p which leads (after careful algebra) 
to an equation of motion for Sf 

.2 



(dt + Vp ■ dx) 6f - np{l + Up] 



m 



2EpT 



dtH = C[5f,p], 



where 



m 



m 



QT2 



T=To 

In this result, we have used the definition of the beta function 



rp2 





Pig) = 



g 



levr^ 



IICa 



-NjTf 



(3.42) 



(3.43) 



(3.44) 



Since the source term in Eq. (j3.42p is proportional to the beta function, the equilibrium distribution 
Up {t,x,p) provides an exact solution of the Boltzmann equation in a conformal theory. In the 
presence of weak conformal breaking, 5f is of order ^ g'^, and Eq. (j3.42p can be solved numerically 
for 6f /g^H using the same procedure as in the previous cases. 

Once 5f is determined, the stress tensor can be found. Here we will follow the analysis of 
Ref. [531] (see also [32l . l54i ]) which determines the appropriate stress tensor in the presence of a 
nontrivial dispersion relation. The stress tensor is 



T^^iX) = -e{TH{X)) + W{Th{X)) - j ^^^SfiX) 



i.e. for i^' 7^ 0, we have 



T^{-e + W) 



To 



(27r)3 Ep 
d^p Sf{oj, k) 



(3.45a) 



H{uj,k). (3.45b) 
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The term in square brackets is r]^^r]ai3G^j^'^^ {oj , k). We note that it is the "tilde" mass that appears 
in Eq. (j3.45ap making the correlator second order in the conformal breaking parameter. Fig. [ll[d) 
shows the bulk spectral function and exhibits both rich hydrodynamic structure and a broad 
response. 

The hydrodynamic prediction is found as follows. Around the equilibrium stress tensor, there 
is a small correction 



T^"'{X) = [e{TH{X)) + V{Th{X))\ U^U^^ + V{TH{X))g^"' {X) + 6T>"' 
which in kinetic theory is given by 



(3.46) 



(3.47) 



However, in the hydrodynamic regime the full stress tensor is parameterized hy e{t,x) and with 



e(t, x) = eiTniX)) + e(t, x) -Co - ^Tcy H{t, x) + e{t, x) , 



U^'iX) = U^iX) + SU^'iX) ^{1- -Hit, x),u'it, x) ) . 



(3.48) 
(3.49) 



Substituting these expressions into the conservation laws (Eq. ()3.16|) and Eq. ()3.17p ) yields the 
linearized equations of motion after careful algebra 



dte + (eo + Vo)dy =^Tc, dtH (l - Sc^) , 
(eo + Vo)dtu' + djT^' =0 , 



(3.50) 
(3.51) 



where the tensor r*-' is 



(3.52) 



In determining these equations we have used the relation, = (co + Vo)/Tcy. In solving these 
equations, we can work to lowest order in the deviation from conformality, — 1/3. Noting that 
C ~ (Cg — 1/3)^ 1361], while the deviations e{t,x) and u^{t,x) are of order (c^ — 1/3), we determine 
T^^{uj, k) to leading order in (c^ — 1/3) for K ^ 



T^^{u,k) 



T^{-e + 3V) 



To 



H{u,k) + {-l + 3cl)e{uj,k) + -iujCH{Lo,k) . (3.53) 



Solving for €{t,x) from Eq. p.50p . substituting this result into T^^, and finally comparing to 
Eq. (|3.45bp . we determine the hydrodynamic prediction for this correlator 



'ntiiy'nai3G'^''"'^{u}, k) = (1- 3c^)^Tc^ 



-cjj'^ — iV^ujk'^ 



Using the thermodynamic result 



d 



— {cgkY + iTsOjk"^ 



d 



(1 - 'iciY Tc. = {3s—-T—] i-e + 3V) 



dT 



(3.54) 



(3.55) 
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the imaginary part can be written 



{U},k) 



2oj 



^ d__ d_\ (c,A;)^r,fc^ 



9C, 



(3.56) 



Thus as /c — )• the first term approaches a delta function 



'n^lur]a|3- 



2oj 



vr 



Cgk) + -5{ijj + Cgk) 



+ % , (3.57) 



in agreement with the previous analysis of Romatschke and Son 49|]. This explains the sharp sound 
pole seen in Fig. 

The free streaming Boltzmann equation does not provide a good description of the asymptotic 
solution at large k and uj. This is because the large k and large uj limits do not commute with the 
p — 7- limit; the bulk channel is sensitive to these soft momenta. Indeed if we neglect the collision 
term in Eq. <^Mi), the "solution", 

"2 (jjHnpil + Up) 



6fiu;,k) 



m 



2EpT oj — Vp ■ k + ie ' 



(3.58) 



does not obey the boundary condition limp^o x(p) = 0, and the T^^ computed with this "solution" 
is infrared divergent. Thus the free result is not shown in Fig. [2|^d). 

IV. COMPARISON WITH SECOND ORDER HYDRODYNAMICS 



In the long wavelength limit the response of the linearized kinetic theory discussed in Section [TI] 
can be described with linearized hydrodynamics, which is a systematic expansion in gradients. 
At leading order in the coupling, the microscopic dynamics described by this kinetic theory is 



conformal, and thus the appropriate hydrodynamic theory is conformal hydrodynamics [50|, l55l ]. 
Beyond leading order in the coupling, there are corrections to the kinetic theory which break scale 
invariance, and a more general non-conformal hydrodynamic theory must be used to characterize 
the long wavelength response of kinetic theory and gluodynamics more generally [56|]. Indeed, in 
Section IIIIEI we determined the subleading non-conformal corrections to kinetic theory due to the 
scale dependence of medium masses and used this result to determine the bulk response function. 
Subsequently we analyzed this response with non-conformal hydrodynamics at first order. 

In this section, we will first limit the discussion to leading order kinetic theory where the 
microscopic dynamics is conformal. We will analyze the sound and tensor channels with conformal 
hydrodynamics through second order following the fundamental work of Baier, Romatschke, Son, 
Starinets, and Stephanov (BRSSS) 50(]. The goal here is to determine the k and oj where the second 
order theory ceases to be a good description of the dynamics. Subsequently we will analyze the 
bulk channel with non-conformal linearized hydrodynamics through second order and determine 
the relevant non-conformal transport coefficients. 



A. The sound mode and conformal second order hydrodynamics 



In conformal linearized hydrodynamics, the dissipative part of the stress tensor, which is con- 
formally invariant and is second order in derivatives, is given by [50(] 

+ K \R^^'''^ - 2uc,R"'^'"'^^U(s] , (4.1) 



vr: 



Wit 



2u„i2°<^^>^M/3 - 2 (V'V InT) 
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FIG. 2: The spectral density p{oj) = — 2IniG_R(aj, fc) for the (a) shear mode G|f^^(w, k), (b) the sound mode 
G^^^{uj,k), (c) the tensor mode G^y^y{uj,k), and (d) the bulk mode rj^y-qapG^"^ {uj,k). The sohd hues 
show the complete results, while the dotted lines show the expectations of the free Boltzmann equation. The 
variables uj and k are measured in units of ^a/T, with = 5^C'A?^il)/87rlog(T/m£)). The shear viscosity 
is v/ie+p) = 0.4613 T/^A so that w = 0.5,1.0,2.0,4.0 corresponds to uj r]/[{e + p)cl] ~ 0.7,1.4,2.8,5.6, as 
chosen in Fig. 21 



where i?"'^'^'^ {R^'^) is the Riemann (Ricci) tensor, and Tjt and k are new transport coefficients. 
Using the zeroth order equations of motion and the conformal dependence of i] on temperature 
■q (xT^, BRSSS rewrite the constituent relation as a dynamical equation for tt'^'^ 



IT' 



TT' 



(4.2) 



We will 



This equation is similar to the phenomenological model of Israel and Stewart |57l . 
compare the static theory, which consists of V^T'^'^ = and a constituent relation (Eq. (14. ip ). to 
the dynamical theory, which consists of V^jT^''' = and a dynamical equation (Eq. (j4.2p ). We 
will see that although the two theories agree in the limit of validity, neither really reproduces the 
structure of kinetic theory for uj,ck > 0.8 [r//(eo + Vo)c^]~^, though the dynamical theory fairs 
better for larger uj. 
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In conformal second order hydrodynamics, the Green function of the tensor channel is deter- 
mined by solving the equations of motion in the external field hxy{t,z). Following the procedure 
described above, we have [50(] 



1 



(4.3) 



When a; = 0, the source term of the Boltzmann equation is zero (see Eq. (13.13P ). while G^^^^(0, k) = 
Therefore, k = in a theory based on the conformal Boltzmann equation to this order 



55l |. At higher orders k is non-zero. Indeed, k is determined by the k dependence of the static 
succeptibility 



d^X e^*'-^^(t) {[T''y{t, x), r^'^(0, 0)])' 



■ — tvk 
2 



(4.4) 



which may be calculated with straightforward perturbation theory. For pure glue this calculation 
has been done in perturbation theory and the result is non-zero, k = dAT^ /18 [491]. Nevertheless, 
we see that k is of order ~ T^, and is significantly smaller than the other second order transport 
coefficient r]Tn ~ T'^/g^ which can be determined by the linearized Boltzmann equation to this 
order. Specifically, we extract riTj^ by examining the real part of the response function in the limit 
w — )• at A: = 0. For Nc = 3 and various numbers of flavors in a leading log approximation we 
have: 








2 


3 


Tn/iv/sT) 


6.32 


6.65 


6.46 



The details of the multicomponent plasma are presented later in Section IV A[ r,r has been deter- 
mined previously in a complete leading order calculation in Ref. [55[. 

Now we will calculate the retarded Green function of the sound channel in two different ways. 
The conservation laws read 

1 



dt€ + (co + Vo)d,u' 



-(eo + Vo)dth, 



(eo + ro)dtu' + c%e + d^7r' 
In the static theory, the constituent relation is 



0. 



2 2 2 



4 cj 



I ^ 3 * 



while in the dynamic theory, Eq. (|4.2p becomes 

2 



rjdthz 



-K,d?hz 



1 

Using the static formulation, we solve the equations of motion and find 



Gr'(w,A:) = (eo + Vo) 



clk"^ + iVgUok'^ — T-jJ^gC^k^ 



In the dynamic theory, we find 

'■{uj,k) = {eo + Vo) 



Gzzzz I 
R 



iT gOJ — iTT^cjuj — o k/ (Co + Vo)u}' 



c^k"^ + iTsiok'^ + iTT^c^Lok'^ 



(dynamic) 



(4.5) 
(4.6) 

(4.7) 

(static) . 

(4.8) 

(4.9) 
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The dispersion relations for the static and dynamic theories are 

- clk^ + iTsOjk'^ - T^TsC^k^ =0 (static) , (4.10) 
u — Cgk +iTsUjk + iT-j^CgUjk — iTt^uj =0 (dynamic). (4-11) 

In the static theory, the dispersion relation has only two solutions 

00 = ±csk - '-Tsk^ ^ Y (i^ - ^-^«) + ^(^') • (^-12) 

By contrast, the dispersion relation in the dynamic theory has the two physical solutions of 
Eq. (|4.12p . and an extra solution 

u = -— + 0{k^). (4.13) 

Since oj remains constant as fc — >• 0, this last solution lies beyond the hydrodynamic approximation 
5C|. 

Fig. [3] compares the full spectral density of sound channel with first and second order hydro- 
dynamics. For u],ck < 0.35 [ri/{eo + Vo)c'i]~^, first order hydrodynamics does a reasonable job 
at capturing the dynamics. The second order theory does a good job for this entire range, and 
continues to work qualitatively until 

u,ck:^ 0.7 [i]/{eo + Vo)cl]~\ (4.14) 

For still larger k, the dynamic second order theory becomes too reactive, while the static second 
order theory becomes too diffusive. Nevertheless, the dynamic theory seems to capture some aspects 
of the high frequency response better than the static theory. All hydrodynamic simulations so far 
have been based on the dynamical theory, which is hyperbolic and causal 0]. 

B. The bulk mode and non-conformal hydro dyanamics at second order 

Now we will perform a similar analysis of the bulk mode, which must be analyzed with non- 



conformal second order hydrodynamics 56|]. In non-conformal linearized hydrodynamics the shear 
strain to second order is 

^^^u ^ ^/.u ^ {Daf"") + K \R^^"'^ - 2uaR''^^"''^^u,3\ + K*2n„i?"<'^'^>'^n'^ , (4.15) 

while the bulk tensor can be written 

n = Hi + CTnD{V ■ u) + + if^u^upW^P . (4.16) 

We will show that the non-conformal couplings to the Rt^'^°'l^ (i.e. k*,^5,^6) all vanish to 
the order considered here. As with k discussed above these couplings are determined by static 
susceptibilities. To see this, we turn on a static gravitational field 

g^^ix) = (1 + H{x))7]^^ + diag(0, h{x), h{x), h{x)) . (4.17) 

Substituting this form into the hydrodynamic equations, we compute a particular combination of 
the stress tensor components 

{2T''{k) - (r^^(A:) + Tyy{k))) = [2K*k'^] H{k) + [Kk^] h{k) . (4.18) 
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FIG. 3: (Color Online) The (a) real and (b) imaginary parts of the retarded Green function in the sound 
channel, G^'^^{uj,k) . The thick solid lines of various colors show the full numerical results from the 
Boltzmann equation; the thin dashed-dotted black lines show the prediction of the first order Navier-Stokes 
equations; the dashed lines show the prediction of the static second order theory (where tt^'' is determined by 
the constituent relation, Eq. (|4.ip ): the dotted lines show the prediction of the dynamic second order theory 
(where n^'^ is determined by a relaxation equation, Eq. (|4.2p '). The shear viscosity is ri/{e+p) = 0.4613 T/ ha 
so that w = 0.1, 0.2, 0.3, 0.5 corresponds to w 7?/[(e + p)cl] ~ 0.14, 0.28, 0.42, 0.7 . 

Similarly, we define 

auik(t, x) = 3ciT%{t, x) + T\{t, x) , 

rjL 

as was motivated in Ref. [53|, and note that in the static graviational field described above 



(auik(fc)) 



H{k)+[6^5k^] h{k). 



(4.19) 



Thus, all of the non-conformal couplings to the curvature tensor are related to the static suscepti- 
bilities 

-i j d^x e'^ '^eit) { [2T''{t, x) - T^'^'xt, x) - Tyy{t, x) , r^(o, o)] )' = 2 K*k'^ (4.20) 

~i J d^Xe'''-^e{t){[O,^,,{t,x),T>^^{0,0)]y = 9^5k' -^Cek' (4.21) 
-i I d^X e^*'-" e{t) (auik(i, x) , T^O, 0))' = 6^5^:' (4.22) 



These may be computed with straightforward perturbation theory as was done for k in Ref. [49(]. 
Further, since = ^^G^yG^^ , every insertion of brings at least two powers of g. Thus we 
estimate that 

z^* = = ^6 = + 0{g^) (4.23) 



In the Boltzmann equation, this must be considered zero to the order we are working. Indeed, we 
see from Eq. (13.13P and Eq. (|3.42p that the sources for 5f induced by H{x) and h{x), 



np{l + Up) 



m 



2EpT 



dtH , and np{l + 



dth, 



(4.24) 
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vanish for time independent gravitational fields. 

To determine the last remaining coefficient rn, we will consder a correlation function of 

Obulk(il 

r]^.Gl^\uj, k) = -il d^X e*-*-*'=-e(t) ( [a.i.(t, x),T>^^{0, 0)] > , (4.25) 

which can be constructed by turning on a gravitational field g^i, = (1 + H{t, x))rj^i, and evaluating 
(auik(t,a3)) 

(auik(w, k)) = -ir7^,G^'"^(u;, A:) . (4.26) 

At = 0, we substitute gfj_u = (1 + H[t))rj^y into the second order non-conformal hydrodynamic 
equations, and determine (C'buik('^, 0)) when ■^s = = 

(auik(t^, 0)) = - ^ [-9<w + 9CTna;2] ^(^^ _ (4^27) 

The quantity in square brackets is the hydrodynamic prediction for the response function and 
should be valid at small w. In kinetic theory we turn on g^^y = (1 + H{t, x))r]^y and measure the 
response C'buik(w,fc) as described in Section Hill 

Obuik(w, /c) = - ^ 

Comparing the functional form of our numerical results from kinetic theory to the hydrodynamic 
form at = 0, we determine rn 



Nf 





2 


3 




0.510 


0.548 


0.554 



(4.29) 



Thus we see that the relaxation time of bulk perturbations is similar to the relaxation time of shear 
perturbations. 

V. EXTENSION TO QUARKS AND SPECTRAL DENSITIES OF 

So far we have only discussed and simulated the pure glue theory. In QCD, the quarks carry 
nearly half of the entropy. For this reason, it is important to extend the analysis to include quarks. 
However, we will see that to ~ 10% accuracy, the overall shape of the T^'^ spectral functions given 
above is unchanged. After including quarks, we determine the current-current correlation function 
in high temperature QCD. 

A. Extension to quarks 

When quarks are added to the mix, the leading log Boltzmann equation becomes somewhat 
more involved. Each species is expanded as follows 

<5r = np(l±np)x"(p). (5.1) 



d'^p 3cfm^ - (1 - 3cf )p^ 5f{u, k) 



E„ 



H{uj,k)/2 



H{uj, k) 



(4.28) 



21 



The collision operator is best expressed in terms of the sum of the fermion and the anti-fermion 
distribution functions, Sf^^^ = (5/^ + 5f^, and the corresponding difference, 5f'^~^ = Sf^ — 6f^. 
The gluon distribution evolves according to 



r 9 _|_ r9 _|_ rs 

'-'FPloss ^ '-'FPgain ^ ^qg 



(5.2) 



{dt + vp-d^)6f9{t,x,p) 

Similarly, the q + q distribution function obeys the equation 

{dt + Vp ■ d^)5r^Ht, x,p)= [(C^p,_ + 4p,_) + (C^pg,,„ + 4pg,, J + (q, + C|,)] . (5.3) 

The corresponding Boltzmann equation for the fermion difference is simpler since the gain term 
for the Fokker Planck evolution cancels in the difference 



The ingredients here are the following: 

1. The Fokker-Planck evolution loss term is 

d 



q 

FPloss 



^FPloss) + i^qg ^Ig)] 



(5.4) 



d 

Cfpioss[x,p] = Tua^ {np{l± np) — 



9p* 



np(l ± Up) 



where the species dependent drag coefficient is 



dt 



and the Debye mass 



9" .2, 



with 



/"a 



Stt 



log 



-) 

mo J 



(5.5) 



(5.6) 



m 



2 _ ^CRa 



<(1±<) g^T^ 



{Nc + NfTp) 



T 3 
Cpta is the quadratic Casimir of species a. 
2. The gain terms are also similar to the previous section. The total diffusion current is 

a " 

which only involves the sum of fermion and anti-fermion flavors. The gain terms are 



(5.7) 



(5.8) 



^FPgain TmldA 



p2 Qp 



P^np{l ± Hp) 



dt ' Tml)dA 



d_ 
dp 



np{l ± Up) 



dP 

~dt 



(5.9) 



where dE/dt and dP /dt axe the total work and momentum transfer per volume on the hard 
particles, i.e. Eqs. (]2.7p and (j2.8p but with the current given by Eq. (j5.8p . 



3. When quarks are present there are additional processes that transform quarks and anti- 
quarks to gluons (and vice versa) by scattering off the fermion mean field. For the fermion 
sum, the collision operator for this process is 



^qg ^ ^qg 



P 



[X%p) + x'^{p)-'2x'{p)\ , 



(5.10) 
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where we have defined for subsequent use 



l = ^^\o^{T/mn), iBF = j^ ^ = (5.11) 

The analogous gluon colhsion operator entering in Eq. (j5.2p is 

f 



where the sum is over the light quark flavors. 

4. For the fermion difference, we note that the Fokker-Planck evolution is the same as before, 
but the gain terms cancel when the difference is taken. The quarks and gluons scattering off 
the fermion mean field (i.e. diagrams D and E in Appendix lA 2\ yield the collision term 

(1 + n^) , _ , 

Kbf P Jk k ^ 
where the first line is the loss term and the second line is the gain term. 

Using these formulas, the spectral functions computed in the last section can also be computed 
in multi-component plasmas - see Appendix [B] for details. We have found that when the response 
functions are expressed in terms of appropriately scaled kinematic variables 

^="(i:TkM' '='^'(i:TkM' ^'-''^ 

the spectral densities are essentially unchanged in all channels. In Fig. Hl^a), we show the bulk 
spectral function in terms of these scaled variables for three flavors and pure glue. The relative 
agreement between these curves indicates the dominance of the Fokker-Planck evolution in deter- 
mining the shape of the response functions. In the next section, we will discuss spectral densities 
of J'*, where a similar scaling is observed if (jJ and k are scaled by the diffusion coefficient, as seen 
in Fig. mb). 



B. Spectral densities of 

For simplicity and definiteness we will compute the current-current correlator of net strangeness. 
Since in a leading log approximation the susceptibilities and correlators are diagonal in fiavor space, 
the flavor and electromagnetic spectral densities are trivially related to this result. To determine 
the strangeness response function, we will probe the Boltzmann equation with a fictitious fiavor 
gauge field that couples only to the strangeness current. 

In the framework of linear response, the average current in the presence of an external field is 

(J^(X))^ = +i j d^y0(xO-yO)^[jM(x),r(y)])^,(y). (5.15) 
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FIG. 4: (a) The bulk spectral function for three flavors compared to the pure glue theory. In this figure, 
?7/(eo + Vo) is 0.917 r///A for Nf = 3 and 0.461 T / ^jla for Nf ^ 0, so that the k values for TV/ = coincide 
with Fig. [5] (b) The longitudinal current-current spectral function for three flavors and the quenched 
approximation. In this figure D = Q.MAT/ ^ip for Nf = 3, while D = 0.852r/^F for Nf = 0, so that the k 
values for Nf = coincide with Fig. [S] The results are similar in the other channels. 



In deriving this result, the definition of the current, J'^(X) = 6S/6A^{X), is used to specify the 
interaction Hamiltonian, Hi^t = — f d?x J^^A^. In Fourier space we define 

G^^(a;, k) = -i j d^X e'^'-'^-'^eit) ([J^(t, x), r{{), 0)]) , (5.16) 

and conclude that 

{r{oj, k))^ = -G^^{u:, k)A,{oj, k) . (5.17) 

Taking k along the z direction, there are two independent correlators, G^j^{uj,k) and G^{uj,k). 

To determine the Boltzmann equation in the presence of an external field, we note that the 
Lorentz force on a charged particle is J-^ = QaF^^v^ , which leads to the Boltzmann equation for 
the strangeness excess 

(p^d^ + QaF^^y,^^ r = C'lf^p] , (5.18) 

where Qs is one for strange quarks, minus one for anti-strange quarks, and zero for all other 
species. Turning on a weak gauge field A^ = (0, A) disturbs the system from equilibrium through 
the linearized Boltzmann equation 

{-iLO + ivp-k)5r{uj,k)-iujnpil-np)QaAi^=r[6f,p]. (5.19) 

We see that the gauge field does not disturb the fermion sum 6f'^^'^, and only disturbs the fermion 
difference 

{-iuj + ivp ■ k)6f-'{uj, k) - iujnp{l - Up) 2Q,Ai^ = C'-'[5f,p] . (5.20) 

Ep 
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0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.5 1 1 .5 2 2.5 3 3.5 4 4.5 



CO = coT/|j,p CO = coT/)j,p 



FIG. 5: (Color Online) The current-current correlator for iVc = 3 and Nj = in the (a) longitudinal 
and (b) transverse channels. Nf ~ corresponds to the quenched approximation. In these k points 
along the z direction, and ^j.f = g^CF'rnjj\og{T/mD)/8Tr is the drag coefficient of a quark in a leading 
log approximation. The diffusion coefficient is _D = 0.852T/ ^f, so ui = 0.5,1.0,2.0,4.0 corresponds to 
LoD/c^ ~ 0.42,0.85,1.7,3.4, as chosen in Fig. The thin dotted lines corresponding by color show the 
results from the free Boltzmann equation. 



The full specification of the collision operator is given in Eq. (I5.13p . and the numerical 

procedure used in the previous section is generalized to solve for 5f'^~^ in Appendix |Bj Once 6f is 
determined (or 6f /QgAz in practice), the current can be determined in a straightforward fashion 

Through Eq. (I5.17p . this determines the current-current response function, which is illustrated in 
Fig.El 

The spectral density for free Boltzmann equation explains the structure of the correlators at 
large lo and k. Following what is by now standard procedure, we find: 

,5,22a) 

and we exhibit these free solutions in Fig^5]with dashed lines. At zero k, the current-current 
correlator has been computed previously [29 1. 

In the long wavelength limit, the current-current correlator is given by the diffusion equation 
which we will develop through second order. The current is expressed in terms of gradients of the 
net strangeness ns{t,x) and the gauge field A'^{t,x). For a linearized theory invariant under parity, 
the current to second order in the derivative expansion must have the following form 

Ji = -D d'us + aE' - (aTj) dtE' + kb (V x . (5.23) 

Here D is the diffusion coefficient, a is the conductivity, and {(ttj) and kb are new transport 
coefficients. In writing this expression, we have neglected the e^^^u^ and ^d^T which would 
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appear in magneto-hydrodynamics [59| (where is not small) or at finite background chemical 
potential (where is not small). Similarly, we have neglected n^u* which is non- linear in the small 
fluctuations of ns{t, x) and u*(t, x). We also have used lower order equations of motion to recognize 
that dtd^n is actually third order in the derivative expansion. 

In fact, the diffusion coefficient and the conductivity are simply related to each other. To see 
this, we first rewrite the constituent relation in terms of the chemical potential, and include one 
higher order term 

J* = —Dxs 5V + <7£'* - (cTTj) dtE^ + Kb (V X By + [ciXs^td^n + other higher order terms] , 

(5.24) 

where Xs is the static susceptibility, 

dn 1 /" 

Xs = ^ = "^Qlvs— J np{l - Up) = Q'ii'.s-^ ■ (5.25) 

Then we note that a perturbation of the form 

fi{X) + Ao{X) = Const (5.26) 

does not drive the system away from equilibrium, i.e. e~^^^"^^^ is constant. Thus all gradients 
in the constituent relation should involve the combination (9*(^ + Aq). This requirement forces a 
relation between the conductivity and the number diffusion coefficient 

XsD = a, (5.27) 

and specifies the coefficient of one higher order term, ci = (Dtj) . 

Now that the constituent relation is specified, the conservation laws dtJ^ + diJ^ = can be 
solved for J^{uj,k) in the presence of a sinusoidal electric field. Through the constituent relation 
(Eq. (I5.23P ) and linear response (Eq. (I5.17P ). this solution completely determines the form of the 
current-current correlator at small momenta 

G''''{u},k) = - iioa + (aTj) u"^ - KBk"^ . (5.28b) 

When oj = 0, the source term of the Boltzmann equation is zero (see Eq. (j5.20p ). while G^^(0, A;) = 
—Ksk'^- Therefore, = in a theory based on the Boltzmann equation. (As in Section HVl this 
transport coefficient may be non-zero at higher order.) In the limit of a; — t- and A; = 0, the real 
part of the Green function gives the value (crj). We tabulate this transport coefhcient here (for 
Nc = 3 and various numbers of flavors in a leading log approximation) 








2 


3 


rj/D 


3.776 


3.756 


3.748 



The imaginary part of Eq. (|5.28ap describes the rapid "dip" seen at small u and k in Fig. O^a). 
The second order theory also makes a definite prediction in Eq. (j5.28ap for the form of the real 
part of the correlator at small u},k. This prediction is shown in Fig. [6l where it is compared to 
the first order theory and the full Boltzmann equation. The second order theory captures some 
aspects of the high frequency response. 
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FIG. 6: (Color Online) The real part of the retarded current current response function for Nc = 3 and 
Nf ~ 0. The thin dotted lines show the predictions of the first order diffusion equation, while the thick 
dashed lines show the prediction of the second order theory, Eq. (|5.28ap . Both are compared to the full 
Boltzmann equation, uj and ck are in units of tJ^ ~ 0.312 /ip/T. Thus fcrj = 0.15, 0.3, 0.6, 1.2 corresponds 
to ujT/ij.f ~ 0.045, 0.09, 0.18, 0.36 in Fig. O i.e. smaller than the first Fig. Rvalue. 



We have written the "static" version of the second order diffusion equation. To the same order 
of accuracy, we can formulate a dynamic second order theory. Using the first order expression 
J] = aE^ — Dd^Hg, we can rewrite Eq. (j5.23p (with kb = 0) as 

Ji = -D d'us + aE' - Tj dtJi . (5.29) 

This is the canonical form of the telegraph equation, which has been extensively studied (sol . 61], 
but we will not develop this connection further. 



VI. CONCLUSIONS 

Fig. [2] shows our principle results for T^'^ correlations in the shear, sound, bulk, and transverse 
tensor channels. These results are for pure glue theory, but in Section IV Al we determined how 
the Boltzmann equation would be modified when quarks are included. Fig. [J] shows how these 
modifications change the bulk and longitudinal current channels for various numbers of colors and 
flavors; the result is similar in the other channels. Fig. Oshows our analogous results for J'* spectral 
densities in the transverse and longitudinal (i.e. the diffusion) channel. 

In each channel (with the exception of the transverse tensor and transverse current), one sees 
a rich hydrodynamic structure at small momentum. We determined the hydrodynamic prediction 
through second order for these channels by solving the conservation laws in an external gravi- 
tational and electromagnetic fields. The electromagnetic response through second order has not 
been determined previously. Fig. [3] analyzes the sound channel in greater detail and compares the 
results to viscous hydrodynamics at first and second order to determine the range of validity of the 
macroscopic theory. (Fig. [6] shows an analogous study of charge diffusion.) Roughly, second order 
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hydrodynamics "works" up to about 



w, ck ^ 0.7 



(6.1) 



For larger momenta hydrodynamics becomes qualitatively wrong, and one sees a transition to free 
streaming quasi-particles (the dotted curves in Fig. [2] and Fig. [5|). However, near the light cone 
uj = k, the momenta need to be truly asymptotic before the sharp structures of the free theory are 
reproduced by the full result. 

Nevertheless, already at fairly modest momenta, a smeared free result describes our full kinetic 
theory result fairly well. Consequently, it will be difficult to distinguish these curves in Euclidean 
space time as their integrals are the same. The question now is, when these smeared spectral shapes 
are combined with a reasonable model of the high frequency continuum, will the full spectral shape 
remain distinctly different from the completely smooth AdS/CFT results? We plan to make a 
complete comparison to the AdS/CFT and the lattice in future work, in an effort to see if quark 
and gluon quasi-particles can be distinguished in the lattice data. 

In addition to providing definite predictions for the spectral densities, we hope that our analysis 
offers a new perspective on the linearized Boltzmann equation. In Section HIl and Appendix (Aj we 
have shown that close to equilibrium, the diffusive motion of the hard particles excess is described by 
a Fokker-Planck equation, which arises from the 2—7-2 scattering graph. The work done during the 
Brownian process shows up as additional gain terms conserving energy and momentum. These gain 
terms have not appeared explicitly before, and they are essential to computing the spectral density 
or to simulating the Boltzmann equation. Ordinarily such a Fokker-Planck equation would conserve 
particle number. However, Bremsstrahlung processes, which are not logarithmically enhanced for 
typical momenta p ~ T, can not be neglected in the limit p — )■ 0. Consequently, the Fokker- 
Planck equation should be solved with an absorptive boundary condition at zero momentum - see 
Section III A II for further discussion. This has been understood previously through an analysis of 



the bulk viscosity [3a], but has not been widely appreciated. In equilibrium (where particle number 
is constant), the particle loss from the temperature scale to the Debye scale through the absorptive 
boundary condition is compensated by the additional gain terms discussed above. A sample of 
evolution of a non-equilibrium distribution ultimately approaching equilibrium is given in Fig. [TJ 
The analysis and numerical work presented here sets the stage for simulating the jet-medium 
response in a leading log approximation. Work is in progress to extend the analysis further, and 
to simulate the jet-medium response in a complete leading order calculation. 
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Appendix A: Leading log analysis of the linearized Boltzmann equation 



This section closely follows Ref. [4g] and determines the leading-log Boltzmann equation for a 
pure glue theory. We will refer to diagrams A-D following the nomenclature of this work. Relative 
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to this work, the gain terms are exphcitly given, and the final form emphasizes the Fokker Planck 
nature of the resulting theory. 



1. Pure glue 

Starting with Eq. (j2.ip . we linearize around the equilibrium distribution writing 

f{t, x,p) = np + np{l + np)x{t, x, p) . (Al) 
The full non-linear collision integral (including a final state symmetry factor) is^ 

C[fM = - I ^|M|2 {27r)^6\Ptot) [fpfkil + + fk') - fp'fk'il + fp){l + fk)] , (A2) 

Jkp'k' ^ 

and this integral is subsequently linearized yielding 

C[f,p] = - I hM\^i27T)^6HPtot)npn,{l + np,)il + nk^)[x{p) + xik)-xip')-x{k')] . (A3) 

Jkp'k' ^ 



In the pure glue theory, the only diagram is t— channel gluon exchange, diagram A of Ref. 46|. 
The linearized integral can be recast as a variational problem 

C[/,p] = -(2vr)3^/M, (A4) 



with 
/ 



[X]^^ [ \M\^{2'K)^6\Ptot)np7i,{l + np,){l + nj,>)[x{p) + xik)-x{p')-xik')Y • (A5) 

-■^0 Jpkp'k' 

We classify the collision integrals as gain and loss terms 

C[f,p] = -(27r)3^A_ (j[;^]i^^^ + , (A6) 

with 

/Mioss = / |M|2 {27T)'6\Ptot)npnk{l + np,){l + nk') -| [x{p) - xip'f , (A7) 

Jpkp'k' J^D 

and 

^Mgain = / |M|2 {2n)^d\Ptot) npnfc(l+npO(l+nfcO ^ [xip) - x{p')] [x{k) - x{k')] ■ (A8) 

Jpkp'k' J^O 

To extract the leading log, we expand in the momentum transfer^, q = p' — P- In a leading log 
approximation, we can neglect the differences between k' and k and the collision integrals read 



/Mioss = np(l + np)^^^0^nk{l + nk)P^ip, k) , (A9) 
-/Mgain = I ^n.p{l+np)^^^^nkil+nk)P\p,k), (AlO) 



^ We will not write the spin and color information explicitly here. The matrix elements are summed over spins and 

colors associated with k,p' , k' and averaged over the spins and colors associated with p. The distribution function 

fp is defined so that the total number of gluons is 2d a /p. 
* We have assumed that p' is close p and that t is small. For identical particles, there is an additional phase space 

region where k' is close p and u is small. This phase space region gives an equal contribution and this factor of 

two is included into the definition of I'-' . 
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where 



(27r)^(5*(P + K-F -K') \M\^ q'q^ . 



(All) 



p'k' 



We will subsequently show that in a leading log approximation P-'{p, k) evaluates to 



B 



(A12) 



where we have defined the leading log coefficient in terms of the parameters ha and S^b given in 
Eq. dial and Eq. (f230D 



Substituting into the loss term yields 



log {T/ mo) . 



-?'[x]loss = ^Tf^A j n.p{l + np) 



dxjp) dxjp) 

Qpi Qpi ' 



(A13) 



(A14) 



where we have used the rotational invariance of nk{l + n^) and the definition ^b- The gain term 
is handled similarly yielding 



-i[x] 



gam 



''2Cb 

TflA 



/ np{l + np)p- 
Jv 



dxjp) 
dp 



+ 



2Cb 



dxjp) 



+ 



Jpk 



npUkil + np){l + Uk) 



/ np(l + Up) 
Jp 

dxjp) dxjk) £ / dxip) \ f dx{k) \ 
dp' dk'J ^' \ dp^ )\ dk' J 



(A15) 



The last line of this equation is in fact zero. To show this, we note that for the rotationally invariant 
nfc(l + nfc) we have 



nfc(l + nfc) 



jj dxjk) jjdxjk) 



dk^ 



Oy dki / k 



Xik) = 0. (A16) 



This result can be used to interchange i and j in k'dx{k)/dy in /[xlgaim yielding our final result 



i[x] 



gam 



2^B 



/ np{l +np)p- 
Jp 



dxjp) 
dp 



+ 



TjlA 
2^B 



np{l + lip) 



dxjp) 



(A17) 



Taking the variation of the gain and loss terms (as prescribed by Eq. ()A6|) ) yields the linearized 
Boltzmann equation given in Section HIl 



a. A leading log integral 

It remains to show that the integral Eq. (jAlip actually equals Eq. ()A12p . To start, we note 
that since the matrix elements are symmetric in p and k, the integral must have the following form 

P^ip, k) = ai {ffP + fc'fc^') + (12 + k'p>^ + a-i5'^ . (A18) 
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To compute the coefficients of the basis we contract with, for instance, p^p^ , p^k^ and 5ij. 
These integrals will be computed in the next paragraph and yield 



p'f^y 

jii 



2^B 



(1 + COS 9pk) = ai(l + cos^ 9pk) + 2a2 cos 6pk + 03 , 

(1 + cos Opjf.) = 2ai cos 9pk + 02(1 + cos^ 9pk) + 03 cos Opf^ , 

(3 — cos 6pk) = 2ai + 2a2 cos 6pk + Sa-^ . 



(Al9a) 
(Al9b) 
(A19c) 



Solving for 01,02,03 yields the result given in Eq. ()A12p . 

To illustrate how the basis integrals are computed, we will compute p^P^fP . First, we use the 
three momentum delta function to perform the k' integral and shift the integral over p' to an 
integral over q, i.e. J 1 ^ J^- Then we rotate our coordinate system so that p points along the 
z axis and k lies in the zx plane; q is measured with respect to this coordinate system, d^q = 

dq (l{cpq)d(t) . (Here and below, we use the shorthand notation Cpq = cos Opq and Spq = suiOpq.) 
In these coordinates, p, k and q = p' — p are 



P = (0, 0, p) , 

k = {kspk, 0, kcpk) , 

q = [qsqp cos 0, qSqp sin 0, qcqp) . 



(A20a) 
(A20b) 
(A20c) 



The energy conservation 5-function can be written 



6{p + k — p — k') 



1 Cpk 



I 



g (1 - Cpfc)^ + Spfc cos 



cos I 



\ 



(1 - Cpfc)^ + Sp^-COS^ 



1/2 



(A21) 



where we have used p' p + q cos 9qp and k' k + k cos O^q. The averaged matrix element (which 
appears in Eq. (jAlip ) in a leading log approximation is 



1 



.g — IQp^k'^Ug 

where the Mandelstam invariants are 



dA ' 



(A22) 



-{P + Kf = 2pkil-cpk), t = -{P'-Pf 



! (1 - Cpk f 

(1 - Cpkf + s 



Ik cos 02 



Thus 



fP\p,k)fP = - 



^|M|2 2t,5{p + k-p' -k')p-qp-q, 



^g9 C% f dq f dcp Spk 2 

— — — ^ / — / — cos : 

SndA J q J 2-kI- Cpk 



'^l^\og{T/mn){l + Cpk) 



(A23) 

(A24) 
(A25) 



The remaining integrals p^P^k^ and are computed similarly yielding the results quoted in 
Eq. (IATqI) . 
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2. Extension to multicomponent plasmas 



We will be quite brief here since our results are to a certain extent simply a reanalysis of Ref. 46 1 
along the lines of the previous section. 

The Boltzmann equation is recast as a variational problem 



(A26) 



Wl 



th^ 



— 16 jpkp'k' 

The collision integral is classified with gain and loss terms as in the previous section. The channel 
exchange diagrams (diagrams A-C) yield 



ah Jpkp'k' 



where the invariant matrix elements are 



M 



(A29) 



Expanding the matrix elements as in the previous section (but keeping track of the species index) 
yields 



lloss _ ^"^D 



\A-C 



167r 



log{T/mD)Y,g^CR^Ua [ n;{l±n;) 



(A30) 



where we have used the definition of the Debye mass 



171% 



(A31) 



The gain terms are handled as in the previous section 



gam 
A-C 



logiT/mp) 
IGirdA 



Y^a'CR^^a I n;{i±n;)p. 



dx^jp) 
dp 



+ 



log (T/mp) 



i:9'CR^^af^nlil±nl)'-^ 



(A32) 



When fermions are included there is also a Compton and annihilation graph. First we will 
handle the annihilation graph. The annihilation graph substituted into Eq. (|A27p . can be written 



^ There is a slight difference between this and the previous section. In the last section, the matrix elements were 
averaged over the spins and colors of the particle a. Here the matrix element is summed over the spins and colors 
of particle a. 
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as sum of a loss and a gain term 

f 

/ pkp 



„ Jpkp'k' 

X [ix'ip) - X'{p)f + ixHk) - X'ik)f] , (A33) 



f 

nx]r = E ^ / l^^'ll' (2^)'5'(^tot)nX(l + n^)(l + n^J 

Jpkn'k' 



I pkp' i 

X [2(x''(p) - x'{p)){xHk) - X'ik))] . (A34) 

The invariant matrix element is 



\Mf,\' ^ Au.CW (-^] ■ (A35) 



Then 

ff 



where the integral 

-12 ^4 



/(p, fc) = I A/p 2^5(pOt) = log(r/mz,) , (A37) 

is easily performed using the parameterization given in the previous section. Then 

nx]T = ^7 E — (^'^(P) - ^'(P))' ' (A38a) 

where we have used the symbols 7 and (,bf given by Eq. (IS.lip . The gain term is handled similarly 
and yields 

= i:^l i^-ik) - X^ik)) <(1±<1 _ . (A38b) 

For the Compton process, we substitute the matrix element 

\Mll\^ ^ -Av.Clg^j (A39) 

into Eq. (|A27p . and read off the loss and gain terms 
ff 

= Y.\I iM^f f (27r)^(5^(Ptot)n^nf (1 + n^)(l - nl){x\p) - x'{p)f , (A40) 

^ ^ Jpkp'k' 

ff 

= E ^ / \M^!\\27r)^6\P,ot)n^n^{l + n^){l - nf ) (A41) 

„ ^ Jpkp'k' 



X [ix^'ip) - X'{p)){x\k) - x^ik))] . (A42) 
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The matrix element is simplified to 

ff 

2^ 



ff 

a -^P 



< (1 + n.^) ^ 



J r 1 gain 



1 7 

2 '^BF 



< (1 + o 



(x"(p)-x^(p)) 



(A43a) 



(A43b) 



Varying according I\x\ with Eqs. ()A30p . ()A32p . ()A38p . and ()A43p . yields the results quoted in 
Section IVAl 



Appendix B: Numerical solution 



1. Pure glue 



For our numerical solutions we introduce a real harmonic basis. 

1 for m = 

Hlm{p) = NimPl\m\{cOSep) X i ^COSTTK^p for m > , 



(Bl) 



-v/2 sin |m|(/)p for m < 



where Nim is the normalization factor 



Im 



2l + l{l-\m\)\ 
47r {l + \m\)\ 



1/2 



(B2) 



and P;|^|(cos0p) is the associated Legendre Polynomial defined without the Condon and Shortly 
phase. We note the equality 



4tt 



(B3) 



For simplicity, we will first discuss the pure glue theory. The LHS (xp^) of Eq. (IXTHD in the 
harmonic basis becomes 



p^LHS = i-iio6u, + ikClP) N{p)xvm - iioNip)nSimip) 



(B4) 



where we have defined 



N{p) = p^np{l ± Up) , 
and recorded the Clebsch Gordan coefficients for k pointing in the z direction 



Cw - ^1+1,1' 



Ni+i^m V 2/ + 1 

In the source term, % is one of the following 

h 



Nim (l-\m\+l\ Ni„, fl + \m\ 



'H = h^x, — , Ky, CA/3{g)T — , 



.H 



(B5) 



(B6) 



(B7) 
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corresponding to the shear, sound, tensor, and bulk modes respectively. Examining Eq. ()3.13p (for 
the first three modes) and Eq. (|3.42p (for the bulk mode), we have the following translations: 



pZpZ 

2EpT 



-VjJ 



pXpV 



N{p)h 



xy 



+iu) 



m 
2EpT 



N{p)H 



where /T'^CAl^ig) is given Eq. (IB37|) . 
Then Eq. (I3l3l) reads^ 



47r 



15 T 



~ 1 

m 



CAfi{g)Tp j ' 



3r 



(B8) 
(B9) 
(BIO) 
(BU) 



iujSw + ikClD)N{p)xi'm - iujN{p)nSim{p) 



jT^^ph -^^(p) 

op op p'^ 



Xlm H : v47r 



+ 




dN{p) 

Cb ' L dp 

dN{p) 2^^, : 
-^^ + -N{p) 
op p 



" dTy 

dP 
' dt 



3r • (B12) 



Im' 



Our goal is to discretize momentum space so that the problem of finding XimiPn) reduces to 
solving a system of linear equations 



(B13) 



To this end, the radial momenta are discretized, pn = 0.5 Ap + nAp with n = . . . M — 1. for 
numerical purposes we define 



w 



Tk 



and set T = 1 from now on. Then the equation of motion becomes 



(B14) 



(_iw,^,;, + ikCl[i)Nip)Fi,„, + Nip)SUp) 



Fir. 



mO 



6 




dN{p) 
dp 



dN{p) 2^^, ; 

-^r^ + -N{p) 

op p 



Cb,p V 3 

where we use a second order difference approximation for the second derivative 

1 



1 d^ 
-iwHfiA dt 
1 dP\ 
i\NHiJ,A dt 



+ 



—N( 
dp dp 



(B15) 



{Apy 



[N{Pn-,l/2) (FiPn+l) - Fipn)) - A^(p„_l/2) (i^(Pn) " i^(Pn-l))] • 



(B16) 



^ The forms of the dE/dt and dP/dt terms in this equation are chosen so that only N{p), which is an analytic 
function, is differenced. Trial and error showed that this has the fastest approach to the continuum. 
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For the energy and momentum terms, we use a midpoint rule 

1 d^\ Ap a Moipn) 



1 dP\ 
iwTinA dt / 1^ 



d_ 

dp' 



(2^)3 Pn -N{Pn)- 




47r X - Ap 



d_ 

dp 



N{pn 



dp ' 
dp 



N{pn)Fim{Pn) 



(BIT) 
(B18) 



The lattice definitions of are defined so that energy and momentum are conserved 



dN{pr, 



dp 



6 



^ 2^ TTTZ^ P n 



dp Pn 



(27r)3 ^' 

and the derivative of the distribution function is 

dN{pn) _ N{Pn+l/2) - N{pn-l/2) 



1 

6 

1 

6 



dp 



Ap 



(B19) 
(B20) 

(B21) 



The boundary conditions of the difference operator in Eq. ()B16P need to be specified. The 
boundary condition discussed in Section [II A 11 x{p)\p=Q = 0, means that we take Foo(p-i) = 

(po), and F/m(p-i) = -F)m(Po) for l>2. In Section UTAH we wrote down 
a first order differential equation at high momentum, Eq. (j2.18p . In the spherical harmonic basis, 
this equation reads 



{-iyNdw + ikCJJf)N{p)Fi,M + N{p)SUp) = -N{p) 



dF, 



Im 



dp 



M7riV(p) 



1 



d^ 



1 



dP 



An 

—N{p) , 

3 V -iwTifJ.A dt ) 



Cs,i5 \-i\NT-L^A dt J ^B,p 

This first order differential equation leads to the update rule for the upper boundary 



(B22) 



FlmiPM) = FirniPAI^l) " Ap (-iw + ikQ™)F,/„ (pA/^ 



^pSlmiPM- 



Ap^-^VA^ 



1 



dE 



-iwTifiA dt 



Ap 




1 



dP 



-iwTifiA dt 



(B23) 



Im' 



We now wish to write the discretized form as the matrix equation, Ax = b. Examining the 
discretized update rules given in Eq. (IBISI) and Eq. (1B16[) . and the boundary condition given in 
Eq. (lB23p . we see that the appropriate vector bnim is 



bnlm = N{pn)Slm{Pn) + Sn,M- 



(Ap) 



)Slm{Pn) • 



(B24) 



We note that the last 6n,M—i piece arises because in Eq. (lB22]) we are specifying the first derivative 
of the distribution function at high momentum. 

In order to solve the system of linear equations, we used BiCGSTAB algorithm which generalizes 
the conjugate gradient algorithm to non-symmetric matrices [64]. In addition to performing the 
multiplication Ax, a typical BiCGSTAB implementation requires A^x. In the present case, A'^x 
involves simply the replacement w — >• — w and k — )• — k in the equations. This is because only 
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the streaming term is the anti-symmetric part of the full matrix. Finally we should specify the 
preconditioner of the conjugate gradient algorithm. Here we simply take the diagonal matrix 
elements when viewed real matrix: 



^precond ^nlm,n'l'm' 



(B25) 



and this seems to provide satisfactory convergence. After solving for Fim, the stress tensor is easily 
found - for example: 

G'-'%uj,k) liA _5T'%u,k) HA ^g26a) 



-lUJ 



cLaT^ +iujhzx cIaT^ ' 



■2\/y|E Pn F2l{Pn) , (B26b) 



n 



(B26c) 



where the overall factor of two is the spin, and the arrow {=^) indicates the limit k = 0,uj ^ 0. 

We record the final expressions for iGR{uj)/uj for the sound channel, the tensor channel, and 
the bulk channel respectively 



6T''iuj,k) flA ^Tf ,Pn L ^ 



(B27d) 

Notice the pleasing similarity with the source terms in Eq. (jBSp . 



2. Multi-component plasmas 



< = ^= E i^aCaJ^npil±np) = ^(^l + ^^Y (B28) 



Now we consider a multicomponent plasma. We introduce a rescaled Debye mass 

2 g,{q + q)/2 

9^Ca 

where we have also rescaled the quadratic Casimir and the number of degrees of freedom 

Ca = ^. i>a = ^. (B29) 

Explicitly we have i>A = 2, and Vq = 2Nf^, and i>(5+5)/2 = ^Nfdp/dA- 
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Then the total work and momentum transfer per volume are 



1 dP\ 
-i\NHiJ,AdA dt J ^ 



a.(q+q)/2 r- — 



(27r)3 ^" dp 
Ap 



dp 



(B30) 



3 ^ (2^)= 



Pn X 



|iV(Pn,.a)^^^ - ^^Nip.^,Sa)F,Up„ 



(B31) 



The lattice versions of the rescaled Debye mass read 

g,(q + q)/2 



a 

g,(q+q)/2 



-^2 _ . A 47r ^ Ap 



dN{pn, Sg) 

dp 

dN{pn,Sa 

dp 



3 ^ (27r) 

The analogous equations of motion for a = g and a = (q + q)/2 are 
=Ca 



H N{pn,Sa) 

Pn 



(B32a) 
(B32b) 
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dp 
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Im 
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-iwli^AdA dt 



+ -N{p, Sg) 
P 



dP 
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mTifiAdA dt ) -vn'H\xa 



To specify the C^^ terms we define 



7= — = 2 7f 



Cp \ £,BP 



CaJ rh'j ' 



"D 



^BF = 7^ , 
lb 



Nbf = pn^p (1 + 



Then the collision terms are 

^2 



p^ 



-iwH^A 



^{9+?)/2 



gg 



^NBFiPn) [2F/r^-'/^(p„)-2Ff„(p„)] 



(B33) 
(B34) 

(B35a) 
(B35b) 



Finally the expressions for the stress tensor remain valid with the appropriate modifications. For 
example, Eq. ()B27eP becomes 



iu{H/2) dAT^ClP{gY 



a.{q + q)/2 



(2^ 



N{pn-, Sg) 



\CAmPn 



9C 



i^A 



dAT^ClP{gf ' 



Fo^obn.) , (B36a) 
(B36b) 
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where the scaled masses are 



_ _|_ - {q + q}/2^ {q + q)/2 (B37) 



In solving the linear system of equations, the transpose is also needed. When multiplying by A^x, 
it must be realized that the matrix implied by Eq. (jB35|) is not symmetric, and the transpose of 
this equation should be used. Alternatively, Eq. (|B35p can be made symmetric by rescaling 
with ^JU^ and changing the formulae of this section appropriately. 



3. Charge diffusion 

In this section we will outline a procedure to solve Eq. (j5.20p numerically. As before we multiply 
by and the left hand side becomes 

p^LHS = {-iiodiv + ikCJ^) N{p)xi'm{p) - iioN{p)ASira{v) , (B38) 

where A is one of 

A = . (B39) 

In this section, N{p) = p^np{l — Up), refers to the fermion distribution and we have dropped the 
s — s label on Xim^- From Eq. (j5.20p . we determine the corresponding sources: 

-iujnp{l - np)2QsA^—— Sf^ = \ —dudmQ , (B40a) 
hipl V o 

-icoiipil - np)2Q,A^^ ^ Sf„, = \ —6iiSmi ■ (B40b) 
Hjpl V o 

For the net strangeness, we define Fim in analogy with the previous section, Fim.{p) ■ 
Eq. (|5.20p becomes 



-iwA ' 



{-isNSlv + ikCli^) N{p)Fi,ra{p) + N{p)SlM 



C 



F 



dp dp p^ 



Flm{p) 

- 2WbAp)FUp) - ^f^V4^NMp) . (B41) 

t.BF,Q \-i\nA dt ) 

where the charge transfer rate is 

^ 2^^Y.l^^B^^Pn)Foo{Pn)- (B42) 



-i\NA dt ^ (27r) 

We choose a lattice definition of ^bf so that the strange charge is exactly conserved 

^BF,Q = 47r J] ^MbM ■ (B43) 
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Finally, from Eq. (j5.2ip and the susceptibility Eq. (j5.25p we determine the longitudinal and trans- 
verse current current correlators (or more precisely iGfi,{uj)/uj) in convenient units 
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